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Abstract. We describe an algorithm computing the exact value of the mean current, its variance, 
and higher order cumulants for stochastic driven systems. The method uses a Rayleigh-Schrodinger 
' perturbation expansion of the generating function of the current, and can be extended to compute 

covariances of multiple currents. As an example of application of the method, we give numerical 
evidence for a simple relation [Eq. ([5])] between the second and the fourth cumulants of the current 
in a symmetric exclusion process. 
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In phenomena like charge transport in nano-electromechanical systems CO, 0, HI 0, H, 0] 
or in life processes like molecular motors or like ion-transport through a single mem- 
brane channel, one easily reaches energy scales as low as a few ksT |0,S]- Therefore, 
the physics and the chemistry of these small systems must talk about fluctuations, not 
only because they are very much present but also because some processes actually make 
use of these fluctuations. Moreover, the experimental output is often in terms of current 
■ cumulants, which should inform us about important features of the system dynamics. 

A central quantity is the fluctuating current j, which gives the time-averaged number 
of particles or quanta that pass trough a given surface, and its large deviation function 
describing the shape of its probability distribution p(j) ~ exp(— TI(j)) for very 
long time intervals T. One is interested in its full shape, as its tail can contain signatures 
of interesting physics. However, the asymptotic estimate of these tails is problematic 
as they involve fluctuations that become rare for T — > °°. We describe a numerical 
scheme to compute exactly (up to numerical round-off 's) the cumulants of these current 
fluctuations in general small Markovian systems. We emphasize that our algorithm is 
exact and can be systematically implemented to produce cumulants of in principle 
arbitrary order. Practical limitations such as memory storage make it however most 
efficient for the statistical mechanics of systems that are not too large. Finally, our 
method is based on general theoretical considerations and its interpretation involves a 
relation between current and traffic (dynamical activity) fluctuations. 

Mathematically, our approach uses a type of Rayleigh-Schrodinger perturbation ex- 
pansion for the largest eigenvalue of a matrix that is obtained as a perturbation of the 
original Markov generator; a full account can be found in [9]. It can be seen as a 



modification and adaption to classical nonequilibrium models, of a technique devel- 
oped within the framework of full counting statistics for quantum transport II 1011 . see 
also [|IL SEE HQ [El]. Our approach is somewhat complementary to a recent ef- 



ficient algorithm for the (approximate) estimate of the large deviation function UM, [170 



EXAMPLE MODEL 

To illustrate our focus we consider the well studied one-dimensional boundary driven 
simple symmetric exclusion process (SSEP), for which several rigorous results are 



available [18, 19]. A state is represented by an array (r] ; ) of N units, which can either 
be empty (rji = 0) or occupied by one particle (T] ; = 1). A transition takes place when 
a particle moves to a neighboring empty site (with rate 1 ) or if it enters/exits from one 
of the end sites, which are in contact with reservoirs at different chemical potentials (a 
at i = 1 and a' at i = N). For example, a transition from a state r\ to a state t, due to 
the entrance of a particle from the left reservoir (rji = — > r\\ = 1 and the rest of the 
array is unchanged) takes place with rate &(rj, <§) = exp(a/2). Observe that we impose 
the physical condition of local detailed balance, meaning that the rates should obey 

"entropy flux" 



*(tj,§) =a(77,£)exp 



(1) 



with some symmetric prefactor a(r],£) = a(^,7]). Here, we have taken a = 1 and the 
irreversible entropy flux from the left particle reservoir is a per entering particle. We 
expect that systematically a net particle current flows through the system, from the side 
with higher chemical potential. As we follow the path or trajectory T](t) over some 
time-interval t E [0,T] we can read the number of particles that enter from the left 
(left time-integrated particle current) and the number of particles that exit to the right. 
These are of course fluctuating currents, as the path (rj (t)) is random with a distribution 
obtained from the Markov dynamics. For our purposes, the information on the dynamics 
is summarized in the generator L, an M x M matrix with elements L(r],£) =k(r],£) for 
7] 7^ £ and L(r\,r\) = — ^^(r/,^). Exactly because of the nonequilibrium condition 
a 7^ a', the matrix L need not even be diagonalizable. For the final algorithm, that 
involves a departure from the more standard Rayleigh-Schrodinger set-up, as we might 
not have an orthogonal basis of eigenvectors. 



PATH-SPACE IDENTITY 

To understand the theoretical point of departure of the method, it is useful to separate 
the time-antisymmetric part from the time-symmetric part in the path-space distribution 
of the Markov process. Under the condition of local detailed balance (0Qh the time- 
antisymmetric part is directly related to the variable entropy flux and the time- symmetric 



part measures the dynamical activity, called traffic ||20LI21|1 in the system. Because of the 
normalization of the path-space distribution these two sectors have related fluctuations, 
as we will make explicit below in ©. Specifically, we consider ensembles of bonds 
B = {rj — ► £ } that all equally contribute to the same time-integrated mesoscopic current 



Jb = Jo dJB(t). In our example, we can take B\ containing all transitions getting one 
particle into the system and coming from the left reservoir, and B2 consisting of all 
transitions in which a particle leaves the system to the right reservoir. We denote by —B 
the ensemble of reversed transitions, giving rise to an instantaneous current dJ b = — 1 • 
We are then interested in the various moments and correlations between the 7g's. These 
can be obtained from the cumulant-generating function 

g(a) = ^ log(exp £a s / s ) (2) 

1 B 

in the steady state of the system. For example the second (partial) derivatives with 
respect to Gb\,Gb 2 give a covariance between two currents. The point is now that the 
exponent in © can be read as an excess entropy flux Y,b&bJb, whose fluctuations are 
the same as that of a dynamical activity in a Markov model with extra driving. To make 
that last point, we imagine an extra driving by modifying the rates to 

L a (r 1 ^)=k(r 1 ^)e^ 

for some antisymmetric function o, which, allowing some abuse of notation, is 
o(t],£) = ±Ob if (77 , ) € ±5. One should check the local detailed balance condition 
(OQ) to see that some extra entropy flux is imposed. Both the original Markov pro- 
cess (with generator L) and the modified one (with generator L a ) have a path-space 
distribution with action, respectively A and A a , and for which 



cq 



with respect to some fixed equilibrium reference dynamics. If we therefore arrange that 
the time- antisymmetric parts of A and A a differ exactly by the excess entropy flux, we 
keep the excess in dynamical activity (time- symmetric parts): 



-A+A a + Y o B J B = I V(rj(t),a) 

B J 



dt 



for a particular function V on the state space, that also depends on the a. That function 
V(rj) essentially measures the difference in escape rates from 77 for the original process 
with respect to the one modified by the a's. The conclusion is a general path-space 
identity 

1 r T 

s(cr) = -log(exp / V{ri{t),a)dt) a (4) 
1 Jo 

where the last expectation is for the modified steady state. The formula © is more ready 
for asymptotic evaluation as T f +°° since V is a multiplication operator. The algorithm 
must then give a systematic expansion in the Cg of the largest eigenvalue (in the sense 

of its real part) of the matrix if = L a + V(; a) = L + M = L + I B £„>i (as)' 1 ' 
where M is the matrix with non-zero elements [e ±cjB — l]k(r),£) only in ensembles 
±5's, where the rates of L a are different from the rates of L. The details of this 



expansion are in [9]. We just briefly mention that this expansion does not have the 
advantage of dealing with symmetric matrices, as in the original Rayleigh-Schrodinger 
method for quantum mechanical operators. The solution of the problem goes via the 
use of resolvents. Interestingly, the final results include the group inverse G of the 
generator L as a main actor 11221, 12311 . This matrix contains all the informations needed to 
compute the dynamical quantities of interest, like the stationary distribution (p | and the 

cumulants. For example, the variance of a current takes the form Cbb = (p|-^b 1-0 ~~ 

(plJ^GJ^Il) (by |1) we mean a vector of l's). 



ILLUSTRATION 

We end with a discussion of some results for the open SSEP introduced above. For the 



SSEP on a ring, cumulants Q} n > of the current have been computed theoretically Il24ll . 
Formulas for the mean current and its variance are also available for the boundar y-d riven 
SSEP, while formulas for cumulants of order n > 3 are known up to order l/N fltL \l $\ . 
Being one-dimensional and finite, there is only one relevant current in the system, which 
we identify with the entrance of a particle from the left reservoir. Our data, for system 
sizes up to N = 14, at equilibrium with half-filling (a = a' = Ofl perfectly agree with the 
theoretical value: = = (2A^) _1 . In the same conditions, for the fourth cumulant 
we observe 

C(4, ^ c(2, » 2 -W (5) 

for N <U, while to order l/N one has = H, QJ]. The relation © is slightly 
different from the one found for the SSEP at equilibrium on a ring, where the asymp- 
totic values of the second and fourth cumulants of the current are related by 



±(£> (2) ) 2 [24]. It turns out that the addition of a term = (£> (2) ) 2 to yields a good 
approximation of for finite N also out of equilibrium, see the Fig. [Eta). Finally, 
we have tested that the relation © is not valid for a = a' ^ (no half-filling) or for 
a = —a' 7^ (a case of half- filling out of equilibrium), see Fig.QJb). 
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The chemical potentials a and a' correspond to reservoirs with particle density p a = 
exp(a/2)/[exp(a/2)+exp(-a/2)) and p b = exp(a'/2)/[exp(a'/2) +exp(-a'/2)) in ifTlfTl . 
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FIGURE 1. (a) Fourth cumulant (C^) of the open SSEP as a function of a, and the theoretical value 
(g( 4 )) by Derrida and coworkers, for N = 14 and a' = 0. The theoretical + (Q^) 2 is also shown, 
(b) (C^ 2 )) 2 and vs a for N = 14, with a' = a and a' = -a (see legend). 
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